Lineage Analysis of CRISPR and shRNA Screens
Sally Claridge
23 August 2018
Analysis of the Broad’s Avana CRISPR (Broad Institute Cancer Dependency Map 2018, Meyers et al. 2017) and the Broad and Dana-Farber Cancer Institite’s Achilles shRNA (MacFarland et al. 2018, Data Science 2018).
The Avana screen produced results using CERES (Meyers et al. 2017) (GitHub), which generates gene dependency scores from sgRNA depletion scores from gene essentiality screens and eliminates bias arising from the effect of copy number variation on Cas9 DNA cleavage. The lower the CERES score, the higher the likelihood that the gene is essential in the associated cell line. Scores are scaled per cell line such that a score of 0 is the median effect of nonessential genes and -1 is the median effect of common core essential genes.
In a previous version of the shRNA screen, DEMETER (GitHub) was used to compute a dependency score for each gene by using the depletion values from each shRNA to infer the effect of target knockdown (on-target) and of expressing a given miRNA seed (off-target) in each cell line. More negative values indicate increased dependency, while more positive values indicate lower dependency. Zero represents the avergae dependency across all cell lines.
In the new data release, DEMETER2 (GitHub repository), developed by McFarland et al. (2018) was used to analyze the Achilles screen. DEMETER2 expands on DEMETER by including parameters for cell-line-specific screen effects and noisy data, correcting for global differences in shRNA levels across cell lines, pooling data acorss cell lines via hierarchical modeling, and using Bayesian inference to compute uncertainty estimates. Now, a score of zero represents no dependency.
All annotation files (copy number, mutation status, and gene expression) (Consortium and Consortium 2015, Barretina et al. 2012) were downloaded from the DepMap Data Portal.
1 Set up
Load libraries.
library(NMF)
library(ggsignif)
library(ggpubr)
library(rowr)
library(data.table)
library(CePa)
library(plyr)
library(tidyverse)
library(magrittr)
library(ggrepel)
library(matrixStats)
library(parallel)
library(kableExtra)
library(gridExtra)
library(broom)
library(glmnet)
library(devtools)
library(reshape2)
library(caret)
library(bsselectR)
library(ComplexHeatmap)
library(circlize)1.1 Functions
1.1.1 Significance thresholding
adj_signif <- function(df) {
# Takes df from compare_means and creates a signif code column for adjusted p-vals
df$p.signif <- ifelse(df$p.signif == "ns", NA, df$p.signif)
df$p.signif.adj <- ifelse(df$p.adj <= 0.0001, "****",
ifelse(df$p.adj <= 0.001, "***",
ifelse(df$p.adj <= 0.01, "**",
ifelse(df$p.adj <= 0.05, "*", NA))))
df$p.short <- formatC(df$p, format = "g", digits = 2)
df$p.adj.short <- formatC(df$p.adj, format = "g", digits = 2)
return(df)
}1.1.2 Make CRISPR grobs
makeCRISPRgrob <- function(g) {
gene <- filter(crispr_data, Hugo_Symbol == g)
# Mutation status
crispr_color <- as.character(gene$Color)
names(crispr_color) <- gene$Mutation_Status
sig <- filter(crispr_signif, Hugo_Symbol == g)
plot_mut <- ggplot(data = gene, aes(x = Mutation_Status, y = Score, color = Mutation_Status)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha = 0.3, size = 0.7, position = position_jitter(w = 0.05)) +
scale_color_manual(values = crispr_color) +
geom_hline(yintercept = 0, linetype = 2, lwd = 0.3) +
theme_light() +
theme(legend.position = "none") +
labs(x = paste0(g, " Mutation Status"), y = "CERES Score",
title = "Mutation Status",
subtitle = paste0("Wilcoxon test:\n- BH-corrected p-value: ", sig$p.adj.short, "\n- Uncorrected p-value: ", sig$p.short))
# Copy number
lo_cn <- range(gene$Copy_Number[!is.na(gene$Copy_Number)])[1]
hi_cn <- range(gene$Copy_Number[!is.na(gene$Copy_Number)])[2]
mid_cn <- (hi_cn - lo_cn) / 2
plot_cn <- ggplot(data = gene, aes(x = Copy_Number, y = Score, color = Mutation_Status)) +
geom_point(size = 0.5, alpha = 0.5) +
geom_smooth(method = "lm", size = 0.5) +
geom_hline(yintercept = 0, linetype = 2, lwd = 0.3) +
scale_color_manual(values = crispr_color) +
stat_cor(method = "pearson", show.legend = FALSE, label.x = c((lo_cn + mid_cn) / 2, (hi_cn - mid_cn) / 2 + mid_cn), label.y = max(gene$Score)) +
theme(legend.position = "none") +
labs(y = "CERES Score", color = "Mutation Status",
x = paste0(g, " Copy number"), title = "Copy Number",
subtitle = "r: Pearson correlation coeffcient")
# Gene expression
lo_ge <- range(gene$RPKM[!is.na(gene$RPKM)])[1]
hi_ge <- range(gene$RPKM[!is.na(gene$RPKM)])[2]
mid_ge <- (hi_ge - lo_ge) / 2
plot_ge <- ggplot(data = gene, aes(x = RPKM, y = Score, color = Mutation_Status)) +
geom_point(size = 0.5, alpha = 0.5) +
geom_smooth(method = "lm", size = 0.5) +
geom_hline(yintercept = 0, linetype = 2, lwd = 0.3) +
scale_color_manual(values = crispr_color) +
stat_cor(method = "pearson", show.legend = FALSE, label.x = c((lo_ge + mid_ge) / 2, (hi_ge - mid_ge) / 2 + mid_ge), label.y = max(gene$Score)) +
theme(legend.position = "none") +
labs(y = "CERES Score", color = "Mutation Status",
x = paste0(g, " Gene Expression (RPKM)"),
title = "Gene Expression",
subtitle = "r: Pearson correlation coeffcient")
# Cell line lineage
plot_lin <- ggplot(data = gene, aes(x = lineage_name, y = Score, color = Mutation_Status)) +
geom_point(alpha = 0.5) +
scale_color_manual(values = crispr_color) +
geom_hline(yintercept = 0, linetype = 2, lwd = 0.3) +
scale_y_continuous(sec.axis = sec_axis(~ .)) +
coord_flip() +
theme(legend.position = "none") +
labs(y = "CERES Score", x = "Cell Line Lineage",
color = "Mutation Status", title = g)
# Arrange plots
plot <- ggarrange(ggarrange(plot_lin, nrow = 1, labels = c("A")),
ggarrange(plot_mut, plot_cn, plot_ge, nrow = 3,
labels = c("B", "C", "D"), heights = c(2, 3, 3)),
font.label = list(size = 30, face = "bold"),
nrow = 1, ncol = 2, widths = c(3, 2))
return(plot)
}1.1.3 Make lineage plot
makeCRISPRlinplot <- function(g) {
gene <- filter(crispr_data, Hugo_Symbol == g)
sig <- filter(crispr_signif_lineage, Hugo_Symbol == g)
gene <- merge(gene, sig, by = c("Hugo_Symbol", "lineage_name"))
gene <- mutate(gene, lineage_name = reorder(lineage_name, p, mean))
crispr_color <- as.character(gene$Color)
names(crispr_color) <- gene$Mutation_Status
score_range <- abs(range(gene$Score)[2] - range(gene$Score)[1])
round_accuracy <- ifelse(score_range <= 2, 0.25,
ifelse(score_range <= 3, 0.5, 1.0))
plot <- ggplot(data = gene, aes(x = lineage_name, y = Score)) +
geom_point(alpha = 0.5, mapping = aes(color = Mutation_Status)) +
scale_color_manual(values = crispr_color) +
coord_cartesian(y = c(min(gene$Score), round_any(x = max(gene$Score), accuracy = round_accuracy, f = ceiling))) +
geom_signif(data = gene, mapping = aes(xmin = lineage_name, xmax = lineage_name, annotations = paste(p.short, "\n", p.adj.short), y_position = max(gene$Score)), manual = TRUE, tip_length = 0, size = 0, textsize = 3) +
theme(legend.position = "top", axis.text.x = element_text(angle = 70, hjust = 1, size = 10)) +
labs(y = "CERES Score",
x = "Lineage",
color = "Mutation Status",
title = paste0(g, ": CERES score by cell line lineage"),
subtitle = "Lineages sorted by increasing p-value; labels indicate unadjusted p-value / BH-corrected p-value.")
return(plot)
}1.1.4 Make individual heatmaps per cluster
Take in a sorted gene list from a larger heatmap and produce smaller, individual heatmaps:
makeHeatmaps <- function(clust_obj, clust_num) {
chm_plot_cols <- colnames(clust_obj)
chm_plot_rows <- rownames(clust_obj)
chm_plot <- Heatmap(clust_obj, name = "CERES Score",
bottom_annotation = HeatmapAnnotation(df = chm_annot),
bottom_annotation_height = unit(5, "in"),
row_title = "Genes", column_title = "Cancer Cell Lines (Broad IDs)",
col = colorRamp2(c(min(chm_scores), 0, max(chm_scores)), c("purple4", "white", "seagreen4")),
na_col = "black",
column_title_gp = gpar(fontsize = 60, fontface = "bold"),
show_column_names = FALSE,
row_title_gp = gpar(fontsize = 60, fontface = "bold"),
show_row_names = FALSE,
row_dend_reorder = TRUE,
column_dend_reorder = TRUE,
column_order = chm_plot_cols,
row_order = chm_plot_rows,
column_dend_height = unit(4, "in"),
row_dend_width = unit(4, "in"),
width = 3)
chm_ge_plot <- Heatmap(chm_ge[row_order(chm_scores_draw)[[clust_num]],], name = "Gene Expression (RPKM)",
col = colorRamp2(c(min(chm_ge[!is.na(chm_ge)]), 0.25, max(chm_ge[!is.na(chm_ge)])), c("dodgerblue4", "white", "firebrick4")),
na_col = "black",
cluster_rows = FALSE, cluster_columns = FALSE,
show_row_names = FALSE, show_column_names = FALSE,
column_title = "Gene Expression",
column_title_gp = gpar(fontsize = 60, fontface = "bold"),
column_order = chm_plot_cols,
row_order = chm_plot_rows,
width = 1)
chm_cn_plot <- Heatmap(chm_cn[row_order(chm_scores_draw)[[clust_num]],], name = "Copy Number",
col = colorRamp2(c(min(chm_cn[!is.na(chm_cn)]), 2, max(chm_cn[!is.na(chm_cn)])), c("dodgerblue4", "white", "firebrick4")),
na_col = "black",
cluster_rows = FALSE, cluster_columns = FALSE,
show_row_names = FALSE, show_column_names = FALSE,
column_title = "Copy Number",
column_title_gp = gpar(fontsize = 60, fontface = "bold"),
column_order = chm_plot_cols,
row_order = chm_plot_rows,
width = 1)
chm_mut_plot <- Heatmap(chm_mut[row_order(chm_scores_draw)[[clust_num]],], name = "Mutation Status",
na_col = "black",
col = c("cyan3", "darkorchid"),
cluster_rows = FALSE, cluster_columns = FALSE,
show_row_names = FALSE, show_column_names = FALSE,
column_title = "Mutation Status",
column_title_gp = gpar(fontsize = 60, fontface = "bold"),
column_order = chm_plot_cols,
row_order = chm_plot_rows,
width = 1)
chm_all_plot <- chm_plot + chm_ge_plot + chm_cn_plot + chm_mut_plot
pdf(paste0("./plots_18Q3/crispr_heatmap_all_c", clust_num, ".pdf"), width = 100, height = 100, paper = "special")
draw(chm_all_plot, heatmap_legend_side = "right", annotation_legend_side = "right")
seekViewport("annotation_Cell_Line_SSMD")
grid.text("Cell Line SSMD", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Cas9_Activity")
grid.text("Cas9 Activity", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Culture_Type")
grid.text("Culture Type", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Primary_Tissue")
grid.text("Primary Tissue", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Primary_Disease")
grid.text("Primary Disease", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Gender")
grid.text("Gender", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Source")
grid.text("Source", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
dev.off()
pdf(paste0("./plots_18Q3/crispr_heatmap_ge_c", clust_num, ".pdf"), width = 100, height = 100, paper = "special")
draw(chm_plot + chm_ge_plot, heatmap_legend_side = "right", annotation_legend_side = "right")
seekViewport("annotation_Cell_Line_SSMD")
grid.text("Cell Line SSMD", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Cas9_Activity")
grid.text("Cas9 Activity", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Culture_Type")
grid.text("Culture Type", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Primary_Tissue")
grid.text("Primary Tissue", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Primary_Disease")
grid.text("Primary Disease", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Gender")
grid.text("Gender", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Source")
grid.text("Source", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
dev.off()
pdf(paste0("./plots_18Q3/crispr_heatmap_mut_c", clust_num, ".pdf"), width = 100, height = 100, paper = "special")
draw(chm_plot + chm_mut_plot, heatmap_legend_side = "right", annotation_legend_side = "right")
seekViewport("annotation_Cell_Line_SSMD")
grid.text("Cell Line SSMD", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Cas9_Activity")
grid.text("Cas9 Activity", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Culture_Type")
grid.text("Culture Type", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Primary_Tissue")
grid.text("Primary Tissue", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Primary_Disease")
grid.text("Primary Disease", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Gender")
grid.text("Gender", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Source")
grid.text("Source", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
dev.off()
pdf(paste0("./plots_18Q3/crispr_heatmap_cn_c", clust_num, ".pdf"), width = 100, height = 100, paper = "special")
draw(chm_plot + chm_cn_plot, heatmap_legend_side = "right", annotation_legend_side = "right")
seekViewport("annotation_Cell_Line_SSMD")
grid.text("Cell Line SSMD", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Cas9_Activity")
grid.text("Cas9 Activity", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Culture_Type")
grid.text("Culture Type", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Primary_Tissue")
grid.text("Primary Tissue", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Primary_Disease")
grid.text("Primary Disease", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Gender")
grid.text("Gender", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Source")
grid.text("Source", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
dev.off()
}1.2 Cell line converter
This comprehensive cancer cell line information curated by Daniel Charytonowicz.
ccl_converter <- read.delim("./data_munging/cell_line_database_FINAL_MASTER.tsv", sep = "\t", row.names = 1, header = TRUE)
colnames(ccl_converter)[colnames(ccl_converter) == "ccle_cell_line_name"] <- "CCLE_Name"
colnames(ccl_converter)[colnames(ccl_converter) == "broad_id"] <- "Broad_ID"1.3 Cell line metadata
ccl_info <- read.delim("./data_munging/DepMap-2018q3-celllines.csv", sep = ",", header = TRUE, na.strings = c("", NA))
ccl_info$Primary.Disease <- gsub("\\\\", "", ccl_info$Primary.Disease)
ccl_info$Primary.Disease <- gsub("Ewings", "Ewing's", ccl_info$Primary.Disease)
crispr_meta <- read.delim("./data_munging/sample_info_18Q3_crispr.csv", sep = ",", header = TRUE, na.strings = c("", NA))
colnames(crispr_meta)[7] <- "CCLE_Name"
shrna_meta <- read.delim("./data_munging/sample_info_18Q3_shrna.csv", sep = ",", header = TRUE, na.strings = c("", NA))
colnames(shrna_meta)[1] <- "CCLE_Name"dan_test <- select(ccl_converter, CCLE_Name, Broad_ID, cell_line_name)
dan_test <- dan_test[-which(dan_test$CCLE_Name == ""), ]
ccl_test <- select(ccl_info, CCLE_Name, Broad_ID)
test <- merge(ccl_test, dan_test, by = "CCLE_Name", all = TRUE, suffixes = c("_DepMap", "_Daniel"))
test$BroadIDs_match <- ifelse(as.character(test$Broad_ID_DepMap) == as.character(test$Broad_ID_Daniel), TRUE, FALSE)1.4 Cancer Gene Census (CGC) gene list
Select genes from the Cancer Gene Census (CGC). The list was pulled from the International Cancer Genome Consortium (ICGC) data portal (Advaced Search > Genes > Curated Gene Set > Cancer Gene Census).
cgc <- data_frame("Hugo_Symbol" = read.delim("./data_munging/gene-ids-for-set-Cancer Gene Census.tsv", header = FALSE, sep = "\t")[, 2])1.5 MAF file
For mutation calling, get paired gene name and cell line fields in a data frame. For this analysis, we don’t care about how many mutations there are per gene or what type of mutations there are, so I didn’t save more information. I took unique gene-cell line combinations since we only cared about mutation presence/absence. Add a Mutation_Status column denoting all entries in MAF files as mutations present in the associated cell lines.
maf_raw <- read.delim("./data_munging/CCLE_DepMap_18q3_maf_20180718.txt.gz", header = TRUE, sep = "\t")
# Remove silent mutations
maf_df <- maf_raw[maf_raw$Variant_Classification != "Silent",]
# Select columns
maf_df <- unique(subset(maf_raw, select = c("Hugo_Symbol", "Tumor_Sample_Barcode", "Broad_ID")))
colnames(maf_df)[2] <- "CCLE_Name"
# Add Mutation_Status column
maf_df$Mutation_Status <- "Mutant"
# MAF summary table
maf_summ <- maf_raw[, c("Variant_Classification", "isDeleterious")] %>% group_by(Variant_Classification, isDeleterious) %>% tally()
maf_summ$percent <- maf_summ$n / sum(maf_summ$n) * 100
knitr::kable(maf_summ, caption = "Distribution of variant classifications in the MAF file") %>% kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))| Variant_Classification | isDeleterious | n | percent |
|---|---|---|---|
| 3’UTR | FALSE | 5 | 0.0004125 |
| 3’UTR | TRUE | 57 | 0.0047020 |
| 5’Flank | FALSE | 12 | 0.0009899 |
| 5’Flank | TRUE | 77 | 0.0063518 |
| 5’UTR | FALSE | 5 | 0.0004125 |
| 5’UTR | TRUE | 30 | 0.0024747 |
| De_novo_Start_OutOfFrame | TRUE | 5179 | 0.4272193 |
| Frame_Shift_Del | TRUE | 39184 | 3.2323152 |
| Frame_Shift_Ins | TRUE | 23445 | 1.9339942 |
| IGR | FALSE | 13 | 0.0010724 |
| IGR | TRUE | 100 | 0.0082491 |
| In_Frame_Del | FALSE | 5412 | 0.4464396 |
| In_Frame_Del | TRUE | 1 | 0.0000825 |
| In_Frame_Ins | FALSE | 1512 | 0.1247259 |
| Intron | FALSE | 20 | 0.0016498 |
| Intron | TRUE | 376 | 0.0310165 |
| Missense_Mutation | FALSE | 666798 | 55.0046277 |
| Nonsense_Mutation | TRUE | 38258 | 3.1559289 |
| Nonstop_Mutation | TRUE | 1014 | 0.0836456 |
| Silent | FALSE | 387651 | 31.9775988 |
| Splice_Site | TRUE | 41381 | 3.4135473 |
| Start_Codon_Del | TRUE | 91 | 0.0075067 |
| Start_Codon_Ins | TRUE | 74 | 0.0061043 |
| Start_Codon_SNP | FALSE | 1272 | 0.1049282 |
| Stop_Codon_Del | FALSE | 2 | 0.0001650 |
| Stop_Codon_Del | TRUE | 171 | 0.0141059 |
| Stop_Codon_Ins | TRUE | 118 | 0.0097339 |
1.6 Copy number file
cn <- read.delim("./data_munging/public_18Q3_gene_cn.csv.gz", sep = ",", check.names = FALSE, header = TRUE)
# Convert log2 ratios (log2(CN/2)) to CN
cn[2:ncol(cn)] <- lapply(cn[2:ncol(cn)], function(x) 2 * (2 ^ x))
# Remove Entrez gene IDs from colnames
colnames(cn) <- gsub(" .*", "", colnames(cn))
colnames(cn)[1] <- "Broad_ID"
# Melt
cn_melt <- melt(data = cn, id.vars = "Broad_ID", measure.vars = colnames(cn[2:ncol(cn)]), variable.name = "Hugo_Symbol", value.name = "Copy_Number")
saveRDS(cn_melt, "./data_munging/rds/cn_melt_18Q3.rds", compress = "xz")cn_melt <- readRDS("./data_munging/rds/cn_melt_18Q3.rds")1.7 Gene expression file
Gene expression data (Reads Per Kilobase of transcript, per Million mapped reads, RPKM).
ge <- read.delim("./../crispr_lineages_giant_files/CCLE_DepMap_18q3_RNAseq_RPKM_20180718.gct.gz", skip = 2, header = TRUE, sep = "\t", check.names = FALSE)
# Edit columns
ge$Name <- NULL
colnames(ge)[1] <- "Hugo_Symbol"
# Melt
ge_melt <- melt(data = ge, id.vars = "Hugo_Symbol", measure.vars = colnames(ge[2:ncol(ge)]), value.name = "RPKM")
## Split variable column
ge_melt <- with(ge_melt, cbind(Hugo_Symbol, colsplit(variable, pattern = " ", names = c("CCLE_Name", "Broad_ID")), RPKM))
## Remove parentheses around Broad IDs
ge_melt$Broad_ID <- gsub("\\(|\\)", "", ge_melt$Broad_ID)
saveRDS(ge_melt, "./../crispr_lineages_giant_files/ge_melt_18Q3.rds", compress = "xz")ge_melt <- readRDS("./../crispr_lineages_giant_files/ge_melt_18Q3.rds")1.8 CRISPR data and annotations
The latest CRISPR CERES score data (18Q3, August 2018) was pulled from the DepMap Data Portal (Broad Institute Cancer Dependency Map 2018, Meyers et al. 2017).
crispr <- read.delim("./data_munging/gene_effect_18Q3.csv.gz", sep = ",", header = TRUE, check.names = FALSE)
# Remove Entrez gene IDs from colnames
colnames(crispr) <- gsub(" .*", "", colnames(crispr))Merge annotation data:
# Melt CRISPR dataset for merging
crispr_melt <- melt(crispr, id.vars = "Broad_ID", measure.vars = colnames(crispr)[2:ncol(crispr)], variable.name = "Hugo_Symbol", value.name = "Score")
# Merge cell line metadata
crispr_melt <- merge(crispr_melt, ccl_info, by = "Broad_ID", all.x = TRUE)
crispr_melt <- merge(crispr_melt, crispr_meta, by = c("CCLE_Name", "Broad_ID"), all.x = TRUE)
# Merge mutation annotations
crispr_muts <- merge(crispr_melt, maf_df, by = c("Hugo_Symbol", "CCLE_Name", "Broad_ID"), all.x = TRUE)
crispr_muts <- crispr_muts %>% mutate(Mutation_Status = if_else(is.na(Mutation_Status), "Wildtype", Mutation_Status))
crispr_muts$Hugo_Symbol <- factor(crispr_muts$Hugo_Symbol)
# Summarize number of mutant and wildtype cell lines
crispr_muts_summ <- crispr_muts %>% group_by(Hugo_Symbol) %>%
summarize(N_Wildtype = sum(Mutation_Status == "Wildtype"),
N_Mutant = sum(Mutation_Status == "Mutant"))
# Merge test results back into full dataset, which restores information lost in the summarization
crispr_data <- merge(crispr_muts_summ, crispr_muts, by = "Hugo_Symbol")
# Add Color column
crispr_data$Color <- ifelse(crispr_data$Mutation_Status == "Wildtype", "cyan3", "darkorchid")
crispr_data$Color <- factor(crispr_data$Color)
# Cell line lineages
crispr_data <- merge(crispr_data, ccl_converter, by = c("CCLE_Name", "Broad_ID"), all.x = TRUE)
levels(crispr_data$lineage_name) <- sort(levels(crispr_data$lineage_name), decreasing = TRUE)
# Copy number
crispr_data <- merge(crispr_data, cn_melt, by = c("Hugo_Symbol", "Broad_ID"), all.x = TRUE)
# Gene expression (RPKM)
ge_filt <- filter(ge_melt, Hugo_Symbol %in% unique(crispr_data$Hugo_Symbol))
crispr_data <- merge(crispr_data, ge_filt, by = c("Hugo_Symbol", "Broad_ID", "CCLE_Name"), all.x = TRUE)
saveRDS(crispr_data, "./data_munging/rds/crispr_data_18Q3.rds", compress = "xz")crispr_data <- readRDS("./data_munging/rds/crispr_data_18Q3.rds")
crispr_ccl <- data.frame("Broad_ID" = crispr_data$Broad_ID)1.9 shRNA data and annotations
The Achilles shRNA DEMETER score data was pulled from the CTD2 Data Portal (Tsherniak et al 2017).
shrna <- read.table("./data_munging/D2_combined_gene_dep_scores.csv.gz", sep = ",", header = TRUE, check.names = FALSE)
colnames(shrna)[1] <- "Hugo_Symbol"
# Remove Entrez gene IDs from gene names
shrna$Hugo_Symbol <- gsub(" .*", "", shrna$Hugo_Symbol)
# Melt shRNA dataset for merging
shrna_melt <- melt(shrna , id.vars = "Hugo_Symbol", measure.vars = colnames(shrna)[2:ncol(shrna)], variable.name = "CCLE_Name", value.name = "Score")
shrna_melt <- drop_na(shrna_melt)Merge annotation data:
# Merge cell line metadata
shrna_melt <- merge(shrna_melt, ccl_info, by = "CCLE_Name", all.x = TRUE)
shrna_melt <- merge(shrna_melt, shrna_meta, by = "CCLE_Name", all.x = TRUE)
# Merge mutation annotations
shrna_muts <- merge(shrna_melt, maf_df, by = c("Hugo_Symbol", "CCLE_Name", "Broad_ID"), all.x = TRUE)
shrna_muts <- shrna_muts %>% mutate(Mutation_Status = if_else(is.na(Mutation_Status), "Wildtype", Mutation_Status))
shrna_muts$Hugo_Symbol <- factor(shrna_muts$Hugo_Symbol)
# Summarize number of mutant and wildtype cell lines
shrna_muts_summ <- shrna_muts %>% group_by(Hugo_Symbol) %>%
summarize(N_Wildtype = sum(Mutation_Status == "Wildtype"),
N_Mutant = sum(Mutation_Status == "Mutant"))
# Merge test results back into full dataset, which restores information lost in the summarization
shrna_data <- merge(shrna_muts_summ, shrna_muts, by = "Hugo_Symbol")
# Add Color column
shrna_data$Color <- ifelse(shrna_data$Mutation_Status == "Wildtype", "cyan3", "darkorchid")
shrna_data$Color <- factor(shrna_data$Color)
# Cell line lineages
shrna_data <- merge(shrna_data, ccl_converter, by = c("CCLE_Name", "Broad_ID"), all.x = TRUE)
levels(shrna_data$lineage_name) <- sort(levels(shrna_data$lineage_name), decreasing = TRUE)
# Copy number
shrna_data <- merge(shrna_data, cn_melt, by = c("Hugo_Symbol", "Broad_ID"), all.x = TRUE)
# Gene expression (RPKM)
ge_filt <- filter(ge_melt, Hugo_Symbol %in% unique(shrna_data$Hugo_Symbol))
shrna_data <- merge(shrna_data, ge_filt, by = c("Hugo_Symbol", "Broad_ID", "CCLE_Name"), all.x = TRUE)
saveRDS(shrna_data, "./../crispr_lineages_giant_files/shrna_data_18Q3.rds", compress = "xz")shrna_data <- readRDS("./../crispr_lineages_giant_files/shrna_data_18Q3.rds")
shrna_ccl <- data.frame("Broad_ID" = shrna_data$Broad_ID)write(union(as.character(distinct(crispr_data[which(is.na(crispr_data$lineage_name)), c("CCLE_Name", "lineage_name")])$CCLE_Name), as.character(distinct(shrna_data[which(is.na(shrna_data$lineage_name)), c("CCLE_Name", "lineage_name")])$CCLE_Name)), file = "./data_munging//na_ccl_in_converter.txt")2 Wilcoxon tests
2.1 CRISPR
2.1.1 Grouped by gene
crispr_signif <- compare_means(Score ~ Mutation_Status, group.by = c("Hugo_Symbol"), data = crispr_data, method = "wilcox.test", p.adjust.method = "BH")
crispr_signif <- adj_signif(crispr_signif)
crispr_signif <- crispr_signif[order(crispr_signif$p),]
saveRDS(crispr_signif, "./data_munging/rds/crispr_signif.rds")crispr_signif <- readRDS("./data_munging/rds/crispr_signif.rds")
knitr::kable(crispr_signif[1:100, c(1, 5, 6, 7, 8, 10)], caption = "CRISPR Screen: Wilcoxon Test Results Comparing Mutant and Wildtype Cell Lines (Benjamini-Hochberg-corrected p-values: * p <= 0.05, ** p <= 0.01, *** p <= 0.001, **** p <= 0.0001)") %>% kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive")) %>% footnote(general = "Wilcoxon test Benjamini-Hochberg-corrected p-values:\n*: p <= 0.05\n**: p <= 0.01\n***: p <= 0.001\n****: p <= 0.0001") %>% scroll_box(width = "900px", height = "450px")| Hugo_Symbol | p | p.adj | p.format | p.signif | p.signif.adj |
|---|---|---|---|---|---|
| KRAS | 0.0000000 | 0.0000000 | < 2e-16 | **** | **** |
| TP53 | 0.0000000 | 0.0000000 | < 2e-16 | **** | **** |
| NRAS | 0.0000000 | 0.0000000 | < 2e-16 | **** | **** |
| BRAF | 0.0000000 | 0.0000000 | 5.0e-15 | **** | **** |
| PIK3CA | 0.0000000 | 0.0000000 | 3.1e-14 | **** | **** |
| PTEN | 0.0000000 | 0.0000000 | 6.6e-14 | **** | **** |
| PIK3R1 | 0.0000094 | 0.0233470 | 9.4e-06 | **** |
|
| GAS2L2 | 0.0000634 | 0.1265210 | 6.3e-05 | **** | NA |
| CTNNB1 | 0.0000658 | 0.1265210 | 6.6e-05 | **** | NA |
| EML5 | 0.0001279 | 0.2113463 | 0.00013 | *** | NA |
| VHL | 0.0001403 | 0.2113463 | 0.00014 | *** | NA |
| ARID1A | 0.0001465 | 0.2113463 | 0.00015 | *** | NA |
| PROM1 | 0.0001756 | 0.2337946 | 0.00018 | *** | NA |
| OR8D1 | 0.0001913 | 0.2364567 | 0.00019 | *** | NA |
| TPR | 0.0003031 | 0.3497722 | 0.00030 | *** | NA |
| TCERG1 | 0.0004767 | 0.5156610 | 0.00048 | *** | NA |
| NFE2L1 | 0.0006240 | 0.6255120 | 0.00062 | *** | NA |
| ATXN2 | 0.0006506 | 0.6255120 | 0.00065 | *** | NA |
| TAOK2 | 0.0007479 | 0.6412831 | 0.00075 | *** | NA |
| FRMD4B | 0.0008177 | 0.6412831 | 0.00082 | *** | NA |
| CAMK2B | 0.0008694 | 0.6412831 | 0.00087 | *** | NA |
| C6orf99 | 0.0009114 | 0.6412831 | 0.00091 | *** | NA |
| KCNIP4 | 0.0009405 | 0.6412831 | 0.00094 | *** | NA |
| TAAR1 | 0.0009576 | 0.6412831 | 0.00096 | *** | NA |
| GMCL1 | 0.0009659 | 0.6412831 | 0.00097 | *** | NA |
| MYH13 | 0.0010170 | 0.6412831 | 0.00102 | ** | NA |
| MMRN2 | 0.0011690 | 0.6412831 | 0.00117 | ** | NA |
| KRT34 | 0.0011703 | 0.6412831 | 0.00117 | ** | NA |
| OR13C4 | 0.0012173 | 0.6412831 | 0.00122 | ** | NA |
| MAT2B | 0.0012266 | 0.6412831 | 0.00123 | ** | NA |
| CCDC115 | 0.0012698 | 0.6412831 | 0.00127 | ** | NA |
| BPHL | 0.0012789 | 0.6412831 | 0.00128 | ** | NA |
| SPHK2 | 0.0013252 | 0.6412831 | 0.00133 | ** | NA |
| SLC12A4 | 0.0013571 | 0.6412831 | 0.00136 | ** | NA |
| HSPB2 | 0.0014495 | 0.6412831 | 0.00145 | ** | NA |
| ING2 | 0.0014811 | 0.6412831 | 0.00148 | ** | NA |
| HSD3B7 | 0.0015616 | 0.6412831 | 0.00156 | ** | NA |
| C11orf80 | 0.0015622 | 0.6412831 | 0.00156 | ** | NA |
| GNAI2 | 0.0015622 | 0.6412831 | 0.00156 | ** | NA |
| PIGW | 0.0015849 | 0.6412831 | 0.00158 | ** | NA |
| PCDHA8 | 0.0016198 | 0.6412831 | 0.00162 | ** | NA |
| ZNF177 | 0.0016563 | 0.6412831 | 0.00166 | ** | NA |
| NCR3LG1 | 0.0016574 | 0.6412831 | 0.00166 | ** | NA |
| EZH2 | 0.0017923 | 0.6412831 | 0.00179 | ** | NA |
| DSEL | 0.0018013 | 0.6412831 | 0.00180 | ** | NA |
| AKR1C1 | 0.0018066 | 0.6412831 | 0.00181 | ** | NA |
| FCGBP | 0.0018105 | 0.6412831 | 0.00181 | ** | NA |
| CPSF1 | 0.0018272 | 0.6412831 | 0.00183 | ** | NA |
| TLX2 | 0.0018522 | 0.6412831 | 0.00185 | ** | NA |
| LRRC37A2 | 0.0018527 | 0.6412831 | 0.00185 | ** | NA |
| BNC2 | 0.0019504 | 0.6618632 | 0.00195 | ** | NA |
| KPNA6 | 0.0021334 | 0.7082436 | 0.00213 | ** | NA |
| TSC2 | 0.0021801 | 0.7082436 | 0.00218 | ** | NA |
| IL20RA | 0.0022216 | 0.7082436 | 0.00222 | ** | NA |
| ST6GALNAC1 | 0.0023064 | 0.7082436 | 0.00231 | ** | NA |
| IL31 | 0.0023228 | 0.7082436 | 0.00232 | ** | NA |
| CMTM4 | 0.0024753 | 0.7082436 | 0.00248 | ** | NA |
| MYO1C | 0.0025576 | 0.7082436 | 0.00256 | ** | NA |
| CLEC4F | 0.0026138 | 0.7082436 | 0.00261 | ** | NA |
| NFE2L2 | 0.0026667 | 0.7082436 | 0.00267 | ** | NA |
| NEDD1 | 0.0028880 | 0.7082436 | 0.00289 | ** | NA |
| ZNF17 | 0.0029504 | 0.7082436 | 0.00295 | ** | NA |
| DIS3L2 | 0.0029826 | 0.7082436 | 0.00298 | ** | NA |
| MSGN1 | 0.0029996 | 0.7082436 | 0.00300 | ** | NA |
| SMCR8 | 0.0030454 | 0.7082436 | 0.00305 | ** | NA |
| ADORA2B | 0.0030699 | 0.7082436 | 0.00307 | ** | NA |
| RBBP9 | 0.0031074 | 0.7082436 | 0.00311 | ** | NA |
| ACSM2B | 0.0031736 | 0.7082436 | 0.00317 | ** | NA |
| ARFGAP1 | 0.0031910 | 0.7082436 | 0.00319 | ** | NA |
| DSPP | 0.0032378 | 0.7082436 | 0.00324 | ** | NA |
| ACBD6 | 0.0032517 | 0.7082436 | 0.00325 | ** | NA |
| SMPD1 | 0.0033308 | 0.7082436 | 0.00333 | ** | NA |
| KRT84 | 0.0033626 | 0.7082436 | 0.00336 | ** | NA |
| OR51T1 | 0.0034115 | 0.7082436 | 0.00341 | ** | NA |
| CD83 | 0.0034583 | 0.7082436 | 0.00346 | ** | NA |
| ETAA1 | 0.0034842 | 0.7082436 | 0.00348 | ** | NA |
| TDRD5 | 0.0035276 | 0.7082436 | 0.00353 | ** | NA |
| HSF4 | 0.0035806 | 0.7082436 | 0.00358 | ** | NA |
| THADA | 0.0035901 | 0.7082436 | 0.00359 | ** | NA |
| FNTA | 0.0036023 | 0.7082436 | 0.00360 | ** | NA |
| ENTPD7 | 0.0036086 | 0.7082436 | 0.00361 | ** | NA |
| ATR | 0.0036293 | 0.7082436 | 0.00363 | ** | NA |
| LY9 | 0.0036616 | 0.7082436 | 0.00366 | ** | NA |
| OR52E8 | 0.0036899 | 0.7082436 | 0.00369 | ** | NA |
| RNF208 | 0.0036899 | 0.7082436 | 0.00369 | ** | NA |
| NPLOC4 | 0.0037050 | 0.7082436 | 0.00370 | ** | NA |
| TMPRSS13 | 0.0037527 | 0.7082436 | 0.00375 | ** | NA |
| SMARCAD1 | 0.0039350 | 0.7082436 | 0.00393 | ** | NA |
| CHRNA1 | 0.0040001 | 0.7082436 | 0.00400 | ** | NA |
| KRT27 | 0.0040103 | 0.7082436 | 0.00401 | ** | NA |
| GUCA1A | 0.0041050 | 0.7082436 | 0.00410 | ** | NA |
| SLC39A3 | 0.0041308 | 0.7082436 | 0.00413 | ** | NA |
| OR5B21 | 0.0042853 | 0.7082436 | 0.00429 | ** | NA |
| TIGD6 | 0.0043171 | 0.7082436 | 0.00432 | ** | NA |
| MRPS5 | 0.0043738 | 0.7082436 | 0.00437 | ** | NA |
| NXPH4 | 0.0044067 | 0.7082436 | 0.00441 | ** | NA |
| AARS | 0.0044532 | 0.7082436 | 0.00445 | ** | NA |
| DACT3 | 0.0044768 | 0.7082436 | 0.00448 | ** | NA |
| SAMM50 | 0.0044860 | 0.7082436 | 0.00449 | ** | NA |
| MTMR2 | 0.0044979 | 0.7082436 | 0.00450 | ** | NA |
| Note: | |||||
|
Wilcoxon test Benjamini-Hochberg-corrected p-values: : p <= 0.05 : p <= 0.01 : p <= 0.001 ****: p <= 0.0001 |
2.1.2 Grouped by gene and lineage
crispr_signif_lineage <- compare_means(Score ~ Mutation_Status, group.by = c("Hugo_Symbol", "lineage_name"), data = crispr_data, method = "wilcox.test", p.adjust.method = "BH")
crispr_signif_lineage <- adj_signif(crispr_signif_lineage)
crispr_signif_lineage <- mutate(crispr_signif_lineage, lineage_name = reorder(lineage_name, p.adj, mean))
saveRDS(crispr_signif_lineage, "./data_munging/rds/crispr_signif_lineage.rds")crispr_signif_lineage <- readRDS("./data_munging/rds/crispr_signif_lineage.rds")2.2 shRNA
2.2.1 Grouped by gene
shrna_signif <- compare_means(Score ~ Mutation_Status, group.by = c("Hugo_Symbol"), data = shrna_data, method = "wilcox.test", p.adjust.method = "BH")
shrna_signif <- adj_signif(shrna_signif)
shrna_signif <- shrna_signif[order(shrna_signif$p),]
saveRDS(shrna_signif, "./data_munging/rds/shrna_signif.rds")shrna_signif <- readRDS("./data_munging/rds/shrna_signif.rds")
knitr::kable(shrna_signif[1:100, c(1, 5, 6, 7, 8, 10)], caption = "shRNA Screen: Wilcoxon Test Results Comparing Mutant and Wildtype Cell Lines (Benjamini-Hochberg-corrected p-values: * p <= 0.05, ** p <= 0.01, *** p <= 0.001, **** p <= 0.0001)") %>% kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive")) %>% footnote(general = "Wilcoxon test Benjamini-Hochberg-corrected p-values:\n*: p <= 0.05\n**: p <= 0.01\n***: p <= 0.001\n****: p <= 0.0001") %>% scroll_box(width = "900px", height = "450px")| Hugo_Symbol | p | p.adj | p.format | p.signif | p.signif.adj |
|---|---|---|---|---|---|
| KRAS | 0.0000000 | 0.0000000 | < 2e-16 | **** | **** |
| TP53 | 0.0000000 | 0.0000000 | < 2e-16 | **** | **** |
| PIK3CA | 0.0000000 | 0.0000000 | < 2e-16 | **** | **** |
| NRAS | 0.0000000 | 0.0000000 | < 2e-16 | **** | **** |
| BRAF | 0.0000000 | 0.0000000 | < 2e-16 | **** | **** |
| CTNNB1 | 0.0000000 | 0.0000000 | 9.0e-12 | **** | **** |
| SMC4 | 0.0000101 | 0.0222880 | 1.0e-05 | **** |
|
| ATR | 0.0000123 | 0.0237737 | 1.2e-05 | **** |
|
| PTEN | 0.0000197 | 0.0338988 | 2.0e-05 | **** |
|
| HRAS | 0.0000848 | 0.1312886 | 8.5e-05 | **** | NA |
| EBNA1BP2 | 0.0001095 | 0.1541798 | 0.00011 | *** | NA |
| TSC1 | 0.0002015 | 0.2599922 | 0.00020 | *** | NA |
| SCNN1D | 0.0002554 | 0.2924992 | 0.00026 | *** | NA |
| FN1 | 0.0002645 | 0.2924992 | 0.00026 | *** | NA |
| XRCC5 | 0.0003202 | 0.3282418 | 0.00032 | *** | NA |
| CDC5L | 0.0003490 | 0.3282418 | 0.00035 | *** | NA |
| TTK | 0.0003604 | 0.3282418 | 0.00036 | *** | NA |
| PIK3R1 | 0.0004016 | 0.3297938 | 0.00040 | *** | NA |
| CKS1B | 0.0004099 | 0.3297938 | 0.00041 | *** | NA |
| ZC3H18 | 0.0004260 | 0.3297938 | 0.00043 | *** | NA |
| DHX9 | 0.0004580 | 0.3376891 | 0.00046 | *** | NA |
| KRTAP5-5 | 0.0005068 | 0.3566407 | 0.00051 | *** | NA |
| C8A | 0.0005485 | 0.3692326 | 0.00055 | *** | NA |
| BOP1 | 0.0006103 | 0.3936929 | 0.00061 | *** | NA |
| QARS | 0.0006713 | 0.4047993 | 0.00067 | *** | NA |
| SRSF1 | 0.0006798 | 0.4047993 | 0.00068 | *** | NA |
| SCAP | 0.0007855 | 0.4439516 | 0.00079 | *** | NA |
| MOS | 0.0008252 | 0.4439516 | 0.00083 | *** | NA |
| ARHGEF15 | 0.0008378 | 0.4439516 | 0.00084 | *** | NA |
| ACP2 | 0.0008603 | 0.4439516 | 0.00086 | *** | NA |
| MORN1 | 0.0009231 | 0.4609898 | 0.00092 | *** | NA |
| POLE2 | 0.0009879 | 0.4640364 | 0.00099 | *** | NA |
| PDE8A | 0.0010152 | 0.4640364 | 0.00102 | ** | NA |
| SPG7 | 0.0010191 | 0.4640364 | 0.00102 | ** | NA |
| SKOR1 | 0.0010556 | 0.4669316 | 0.00106 | ** | NA |
| OR7A5 | 0.0011375 | 0.4797953 | 0.00114 | ** | NA |
| AAMP | 0.0011677 | 0.4797953 | 0.00117 | ** | NA |
| ARID3C | 0.0011791 | 0.4797953 | 0.00118 | ** | NA |
| NBN | 0.0012289 | 0.4797953 | 0.00123 | ** | NA |
| ATP11B | 0.0012396 | 0.4797953 | 0.00124 | ** | NA |
| CYP51A1 | 0.0013130 | 0.4861161 | 0.00131 | ** | NA |
| HSPA1L | 0.0013187 | 0.4861161 | 0.00132 | ** | NA |
| MPP4 | 0.0014354 | 0.5168174 | 0.00144 | ** | NA |
| TPT1 | 0.0014998 | 0.5224054 | 0.00150 | ** | NA |
| NPBWR1 | 0.0015184 | 0.5224054 | 0.00152 | ** | NA |
| FBXO22 | 0.0015556 | 0.5235587 | 0.00156 | ** | NA |
| MED7 | 0.0016389 | 0.5308797 | 0.00164 | ** | NA |
| ELAVL2 | 0.0016614 | 0.5308797 | 0.00166 | ** | NA |
| NXF3 | 0.0017037 | 0.5308797 | 0.00170 | ** | NA |
| ASPM | 0.0017145 | 0.5308797 | 0.00171 | ** | NA |
| GRM2 | 0.0017676 | 0.5365943 | 0.00177 | ** | NA |
| ANAPC4 | 0.0018051 | 0.5374214 | 0.00181 | ** | NA |
| APC | 0.0019759 | 0.5513917 | 0.00198 | ** | NA |
| ZCCHC14 | 0.0020332 | 0.5513917 | 0.00203 | ** | NA |
| SKIV2L2 | 0.0020494 | 0.5513917 | 0.00205 | ** | NA |
| PGBD2 | 0.0020721 | 0.5513917 | 0.00207 | ** | NA |
| SKI | 0.0020893 | 0.5513917 | 0.00209 | ** | NA |
| EMX2 | 0.0021436 | 0.5513917 | 0.00214 | ** | NA |
| DYNC1I2 | 0.0021988 | 0.5513917 | 0.00220 | ** | NA |
| HIST1H2BO | 0.0022369 | 0.5513917 | 0.00224 | ** | NA |
| TRMT61B | 0.0022560 | 0.5513917 | 0.00226 | ** | NA |
| ACVR1C | 0.0022657 | 0.5513917 | 0.00227 | ** | NA |
| RAB3D | 0.0022689 | 0.5513917 | 0.00227 | ** | NA |
| SOS1 | 0.0022794 | 0.5513917 | 0.00228 | ** | NA |
| TBX3 | 0.0024515 | 0.5791092 | 0.00245 | ** | NA |
| DUSP9 | 0.0024688 | 0.5791092 | 0.00247 | ** | NA |
| LDLR | 0.0026312 | 0.5980596 | 0.00263 | ** | NA |
| RASGEF1A | 0.0026400 | 0.5980596 | 0.00264 | ** | NA |
| MAST3 | 0.0026956 | 0.5980596 | 0.00270 | ** | NA |
| CD74 | 0.0027041 | 0.5980596 | 0.00270 | ** | NA |
| SH2D6 | 0.0029572 | 0.6388422 | 0.00296 | ** | NA |
| EPG5 | 0.0029710 | 0.6388422 | 0.00297 | ** | NA |
| PDHA1 | 0.0031875 | 0.6603827 | 0.00319 | ** | NA |
| COL1A1 | 0.0032063 | 0.6603827 | 0.00321 | ** | NA |
| ATG9A | 0.0032157 | 0.6603827 | 0.00322 | ** | NA |
| NDUFS3 | 0.0032603 | 0.6603827 | 0.00326 | ** | NA |
| SF3B2 | 0.0033202 | 0.6603827 | 0.00332 | ** | NA |
| TNFSF8 | 0.0033271 | 0.6603827 | 0.00333 | ** | NA |
| PDLIM5 | 0.0035183 | 0.6841126 | 0.00352 | ** | NA |
| CTSB | 0.0035350 | 0.6841126 | 0.00354 | ** | NA |
| DLD | 0.0035914 | 0.6864390 | 0.00359 | ** | NA |
| ANKH | 0.0036671 | 0.6923699 | 0.00367 | ** | NA |
| TAF3 | 0.0039064 | 0.7083076 | 0.00391 | ** | NA |
| PLP1 | 0.0039093 | 0.7083076 | 0.00391 | ** | NA |
| TSC2 | 0.0039381 | 0.7083076 | 0.00394 | ** | NA |
| MAPK15 | 0.0039958 | 0.7083076 | 0.00400 | ** | NA |
| EZH2 | 0.0040159 | 0.7083076 | 0.00402 | ** | NA |
| RPL3L | 0.0040861 | 0.7083076 | 0.00409 | ** | NA |
| MCUR1 | 0.0041007 | 0.7083076 | 0.00410 | ** | NA |
| COL9A3 | 0.0041739 | 0.7083076 | 0.00417 | ** | NA |
| ROS1 | 0.0041829 | 0.7083076 | 0.00418 | ** | NA |
| EHBP1 | 0.0042584 | 0.7083076 | 0.00426 | ** | NA |
| URI1 | 0.0043305 | 0.7083076 | 0.00433 | ** | NA |
| WRN | 0.0044242 | 0.7083076 | 0.00442 | ** | NA |
| BMP6 | 0.0044278 | 0.7083076 | 0.00443 | ** | NA |
| PRPF39 | 0.0044681 | 0.7083076 | 0.00447 | ** | NA |
| SIGIRR | 0.0044808 | 0.7083076 | 0.00448 | ** | NA |
| TAS2R9 | 0.0044835 | 0.7083076 | 0.00448 | ** | NA |
| STXBP2 | 0.0046239 | 0.7102940 | 0.00462 | ** | NA |
| SETD1B | 0.0046294 | 0.7102940 | 0.00463 | ** | NA |
| Note: | |||||
|
Wilcoxon test Benjamini-Hochberg-corrected p-values: : p <= 0.05 : p <= 0.01 : p <= 0.001 ****: p <= 0.0001 |
2.2.2 Grouped by gene and lineage
shrna_signif_lineage <- compare_means(Score ~ Mutation_Status, group.by = c("Hugo_Symbol", "lineage_name"), data = shrna_data, method = "wilcox.test", p.adjust.method = "BH")
shrna_signif_lineage <- adj_signif(shrna_signif_lineage)
shrna_signif_lineage <- mutate(shrna_signif_lineage, lineage_name = reorder(lineage_name, p.adj, mean))
saveRDS(shrna_signif_lineage, "./../crispr_lineages_giant_files/shrna_signif_lineage.rds")shrna_signif_lineage <- readRDS("./../crispr_lineages_giant_files/shrna_signif_lineage.rds")3 Heatmaps
3.1 CRISPR
3.1.1 Data management
Make and save heatmap data matrices:
# Set up heatmap data
chm_scores <- t(crispr)
colnames(chm_scores) <- chm_scores[1, ]
chm_scores <- chm_scores[-1, ]
chm_rows <- rownames(chm_scores)
chm_cols <- colnames(chm_scores)
chm_scores <- apply(chm_scores, 2, as.double)
chm_scores <- apply(chm_scores, 2, function(x) (x - mean(x)) / var(x))
colnames(chm_scores) <- chm_cols
rownames(chm_scores) <- chm_rows
saveRDS(chm_scores, "./data_munging/rds/crispr_heatmap_scores.rds")
# Make full grid for genomic features
crispr_grid_ccls <- unique(crispr_data$Broad_ID)
crispr_grid_genes <- unique(crispr_data$Hugo_Symbol)
crispr_grid <- expand.grid("Broad_ID" = crispr_grid_ccls, "Hugo_Symbol" = crispr_grid_genes)
# Gene expression
ge_filt <- filter(ge_melt, Hugo_Symbol %in% unique(crispr_data$Hugo_Symbol))[, c("Hugo_Symbol", "Broad_ID", "RPKM")]
chm_ge_grid <- merge(crispr_grid, ge_filt, by = c("Broad_ID", "Hugo_Symbol"), all.x = TRUE)
chm_ge <- reshape(chm_ge_grid, idvar = "Hugo_Symbol", timevar = "Broad_ID", direction = "wide")
chm_ge_rows <- chm_ge$Hugo_Symbol
chm_ge$Hugo_Symbol <- NULL
chm_ge_cols <- colnames(chm_ge) %>% gsub("RPKM.", "", .)
chm_ge <- matrix(as.numeric(unlist(chm_ge)), nrow = nrow(chm_ge))
rownames(chm_ge) <- chm_ge_rows
colnames(chm_ge) <- chm_ge_cols
saveRDS(chm_ge, "./data_munging/rds/crispr_heatmap_ge.rds")
# Mutation status
chm_mut_grid <- merge(crispr_grid, maf_df[, c("Broad_ID", "Hugo_Symbol", "Mutation_Status")], by = c("Broad_ID", "Hugo_Symbol"), all.x = TRUE)
chm_mut_grid$Mutation_Status <- ifelse(is.na(chm_mut_grid$Mutation_Status), 0, 1)
chm_mut <- reshape(chm_mut_grid, idvar = "Hugo_Symbol", timevar = "Broad_ID", direction = "wide")
chm_mut_rows <- chm_mut$Hugo_Symbol
chm_mut$Hugo_Symbol <- NULL
chm_mut_cols <- colnames(chm_mut) %>% gsub("Mutation_Status.", "", .)
chm_mut <- matrix(as.numeric(unlist(chm_mut)), nrow = nrow(chm_mut))
rownames(chm_mut) <- chm_mut_rows
colnames(chm_mut) <- chm_mut_cols
saveRDS(chm_mut, "./data_munging/rds/crispr_heatmap_mut.rds")
# Copy number
cn_filt <- filter(cn_melt, Hugo_Symbol %in% unique(crispr_data$Hugo_Symbol))
chm_cn_grid <- merge(crispr_grid, cn_melt, by = c("Broad_ID", "Hugo_Symbol"), all.x = TRUE)
chm_cn <- reshape(chm_cn_grid, idvar = "Hugo_Symbol", timevar = "Broad_ID", direction = "wide")
chm_cn_rows <- chm_cn$Hugo_Symbol
chm_cn$Hugo_Symbol <- NULL
chm_cn_cols <- colnames(chm_cn) %>% gsub("Copy_Number.", "", .)
chm_cn <- matrix(as.numeric(unlist(chm_cn)), nrow = nrow(chm_cn))
rownames(chm_cn) <- chm_cn_rows
colnames(chm_cn) <- chm_cn_cols
saveRDS(chm_cn, "./data_munging/rds/crispr_heatmap_cn.rds")Load heatmap matrices:
chm_scores <- readRDS("./data_munging/rds/crispr_heatmap_scores.rds")
chm_ge <- readRDS("./data_munging/rds/crispr_heatmap_ge.rds")
chm_mut <- readRDS("./data_munging/rds/crispr_heatmap_mut.rds")
chm_cn <- readRDS("./data_munging/rds/crispr_heatmap_cn.rds")
# Metadata anotation dataframe
chm_annot <- merge(merge(data.frame("Broad_ID" = crispr$Broad_ID), crispr_meta[, c("Broad_ID", "cell_line_SSMD", "cas9_activity", "culture_type", "primary_tissue")], by = "Broad_ID", all.x = TRUE), ccl_info[, c("Broad_ID", "Primary.Disease", "Gender", "Source")], by = "Broad_ID", all.x = TRUE)
rownames(chm_annot) <- chm_annot$Broad_ID
chm_annot$Broad_ID <- NULL
colnames(chm_annot) <- c("Cell_Line_SSMD", "Cas9_Activity", "Culture_Type", "Primary_Tissue", "Primary_Disease", "Gender", "Source")
chm_annot$Cell_Line_SSMD <- as.numeric(chm_annot$Cell_Line_SSMD)
chm_annot$Cas9_Activity <- as.numeric(chm_annot$Cas9_Activity)3.1.2 K-means clustered
3.1.2.1 Euclidean distance
Code for k-means clustered heatmap:
chm_scores_plot <- Heatmap(chm_scores, name = "CERES Score",
bottom_annotation = HeatmapAnnotation(df = chm_annot),
bottom_annotation_height = unit(5, "in"),
row_title = "Genes", column_title = "Cancer Cell Lines (Broad IDs)",
col = colorRamp2(c(min(chm_scores), 0, max(chm_scores)), c("purple4", "white", "seagreen4")),
na_col = "black",
column_title_gp = gpar(fontsize = 60, fontface = "bold"),
show_column_names = FALSE,
row_title_gp = gpar(fontsize = 60, fontface = "bold"),
show_row_names = FALSE,
row_dend_reorder = TRUE,
column_dend_reorder = TRUE,
clustering_distance_rows = "euclidean",
clustering_distance_columns = "euclidean",
clustering_method_rows = "complete",
clustering_method_columns = "complete",
column_dend_height = unit(4, "in"),
row_dend_width = unit(4, "in"),
width = 3,
km = 5, gap = unit(0.5, "in"))
chm_scores_draw <- draw(chm_scores_plot, heatmap_legend_side = "right", annotation_legend_side = "right")
pdf("./plots_18Q3/crispr_heatmap_euclidean_split5.pdf", width = 100, height = 100, paper = "special")
chm_scores_draw
seekViewport("annotation_Cell_Line_SSMD")
grid.text("Cell Line SSMD", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Cas9_Activity")
grid.text("Cas9 Activity", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Culture_Type")
grid.text("Culture Type", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Primary_Tissue")
grid.text("Primary Tissue", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Primary_Disease")
grid.text("Primary Disease", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Gender")
grid.text("Gender", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Source")
grid.text("Source", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
dev.off()
chm_scores_c1 <- chm_scores[row_order(chm_scores_draw)[[1]],]
chm_scores_c2 <- chm_scores[row_order(chm_scores_draw)[[2]],]
chm_scores_c3 <- chm_scores[row_order(chm_scores_draw)[[3]],]
chm_scores_c4 <- chm_scores[row_order(chm_scores_draw)[[4]],]
chm_scores_c5 <- chm_scores[row_order(chm_scores_draw)[[5]],]
makeHeatmaps(chm_scores_c1, 1)
makeHeatmaps(chm_scores_c2, 2)
makeHeatmaps(chm_scores_c3, 3)
makeHeatmaps(chm_scores_c4, 4)
makeHeatmaps(chm_scores_c5, 5)3.1.3 CGC genes
Code for heatmap of CGC-filtered genes only
3.1.3.1 Euclidean distance
chm_scores_cgc <- subset(chm_scores, rownames(chm_scores) %in% cgc$Hugo_Symbol)
chm_ge_cgc <- subset(chm_ge, rownames(chm_ge) %in% cgc$Hugo_Symbol)
chm_mut_cgc <- subset(chm_mut, rownames(chm_mut) %in% cgc$Hugo_Symbol)
chm_cn_cgc <- subset(chm_cn, rownames(chm_cn) %in% cgc$Hugo_Symbol)
chm_plot_cgc <- Heatmap(chm_scores_cgc, name = "CERES Score",
bottom_annotation = HeatmapAnnotation(df = chm_annot),
bottom_annotation_height = unit(5, "in"),
row_title = "Genes", column_title = "Cancer Cell Lines (Broad IDs)",
col = colorRamp2(c(min(chm_scores_cgc), 0, max(chm_scores_cgc)), c("purple4", "white", "seagreen4")),
na_col = "black",
row_title_gp = gpar(fontsize = 60, fontface = "bold"),
row_names_gp = gpar(fontsize = 10),
column_title_gp = gpar(fontsize = 60, fontface = "bold"),
column_names_gp = gpar(fontsize = 8),
row_dend_reorder = TRUE,
column_dend_reorder = TRUE,
clustering_distance_rows = "euclidean",
clustering_distance_columns = "euclidean",
clustering_method_rows = "complete",
clustering_method_columns = "complete",
column_dend_height = unit(4, "in"),
row_dend_width = unit(4, "in"),
width = 3)
chm_plot_cgc_cols <- colnames(chm_scores_cgc)[unlist(column_order(chm_plot_cgc))]
chm_plot_cgc_rows <- rownames(chm_scores_cgc)[unlist(row_order(chm_plot_cgc))]
chm_ge_plot_cgc <- Heatmap(chm_ge_cgc, name = "Gene Expression (RPKM)",
col = colorRamp2(c(min(chm_ge_cgc[!is.na(chm_ge_cgc)]), 0.25, max(chm_ge_cgc[!is.na(chm_ge_cgc)])), c("dodgerblue4", "white", "firebrick4")),
na_col = "black",
cluster_rows = FALSE, cluster_columns = FALSE,
show_row_names = FALSE, show_column_names = FALSE,
column_title = "Gene Expression",
column_title_gp = gpar(fontsize = 60, fontface = "bold"),
column_order = chm_plot_cgc_cols,
row_order = chm_plot_cgc_rows,
width = 1)
chm_cn_plot_cgc <- Heatmap(chm_cn_cgc, name = "Copy Number",
col = colorRamp2(c(min(chm_cn_cgc[!is.na(chm_cn_cgc)]), 2, max(chm_cn_cgc[!is.na(chm_cn_cgc)])), c("dodgerblue4", "white", "firebrick4")),
na_col = "black",
cluster_rows = FALSE, cluster_columns = FALSE,
show_row_names = FALSE, show_column_names = FALSE,
column_title = "Copy Number",
column_title_gp = gpar(fontsize = 60, fontface = "bold"),
column_order = chm_plot_cgc_cols,
row_order = chm_plot_cgc_rows,
width = 1)
chm_mut_plot_cgc <- Heatmap(chm_mut_cgc, name = "Mutation Status",
na_col = "black",
col = c("cyan3", "darkorchid"),
cluster_rows = FALSE, cluster_columns = FALSE,
show_row_names = FALSE, show_column_names = FALSE,
column_title = "Mutation Status",
column_title_gp = gpar(fontsize = 60, fontface = "bold"),
column_order = chm_plot_cgc_cols,
row_order = chm_plot_cgc_rows,
width = 1)
chm_all_plot_cgc <- chm_plot_cgc + chm_ge_plot_cgc + chm_cn_plot_cgc + chm_mut_plot_cgc
pdf("./plots_18Q3/crispr_heatmap_euclidean_all_cgc.pdf", width = 100, height = 100, paper = "special")
draw(chm_all_plot_cgc, heatmap_legend_side = "right", annotation_legend_side = "right")
seekViewport("annotation_Cell_Line_SSMD")
grid.text("Cell Line SSMD", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Cas9_Activity")
grid.text("Cas9 Activity", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Culture_Type")
grid.text("Culture Type", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Primary_Tissue")
grid.text("Primary Tissue", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Primary_Disease")
grid.text("Primary Disease", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Gender")
grid.text("Gender", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Source")
grid.text("Source", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
dev.off()3.1.3.2 Pearson correlation
chm_scores_cgc <- subset(chm_scores, rownames(chm_scores) %in% cgc$Hugo_Symbol)
chm_ge_cgc <- subset(chm_ge, rownames(chm_ge) %in% cgc$Hugo_Symbol)
chm_mut_cgc <- subset(chm_mut, rownames(chm_mut) %in% cgc$Hugo_Symbol)
chm_cn_cgc <- subset(chm_cn, rownames(chm_cn) %in% cgc$Hugo_Symbol)
chm_plot_cgc <- Heatmap(chm_scores_cgc, name = "CERES Score",
bottom_annotation = HeatmapAnnotation(df = chm_annot),
bottom_annotation_height = unit(5, "in"),
row_title = "Genes", column_title = "Cancer Cell Lines (Broad IDs)",
col = colorRamp2(c(min(chm_scores_cgc), 0, max(chm_scores_cgc)), c("purple4", "white", "seagreen4")),
na_col = "black",
row_title_gp = gpar(fontsize = 60, fontface = "bold"),
row_names_gp = gpar(fontsize = 10),
column_title_gp = gpar(fontsize = 60, fontface = "bold"),
column_names_gp = gpar(fontsize = 8),
row_dend_reorder = TRUE,
column_dend_reorder = TRUE,
clustering_distance_rows = "pearson",
clustering_distance_columns = "pearson",
clustering_method_rows = "complete",
clustering_method_columns = "complete",
column_dend_height = unit(4, "in"),
row_dend_width = unit(4, "in"),
width = 3)
chm_plot_cgc_cols <- colnames(chm_scores_cgc)[unlist(column_order(chm_plot_cgc))]
chm_plot_cgc_rows <- rownames(chm_scores_cgc)[unlist(row_order(chm_plot_cgc))]
chm_ge_plot_cgc <- Heatmap(chm_ge_cgc, name = "Gene Expression (RPKM)",
col = colorRamp2(c(min(chm_ge_cgc[!is.na(chm_ge_cgc)]), 0.25, max(chm_ge_cgc[!is.na(chm_ge_cgc)])), c("dodgerblue4", "white", "firebrick4")),
na_col = "black",
cluster_rows = FALSE, cluster_columns = FALSE,
show_row_names = FALSE, show_column_names = FALSE,
column_title = "Gene Expression",
column_title_gp = gpar(fontsize = 60, fontface = "bold"),
column_order = chm_plot_cgc_cols,
row_order = chm_plot_cgc_rows,
width = 1)
chm_cn_plot_cgc <- Heatmap(chm_cn_cgc, name = "Copy Number",
col = colorRamp2(c(min(chm_cn_cgc[!is.na(chm_cn_cgc)]), 2, max(chm_cn_cgc[!is.na(chm_cn_cgc)])), c("dodgerblue4", "white", "firebrick4")),
na_col = "black",
cluster_rows = FALSE, cluster_columns = FALSE,
show_row_names = FALSE, show_column_names = FALSE,
column_title = "Copy Number",
column_title_gp = gpar(fontsize = 60, fontface = "bold"),
column_order = chm_plot_cgc_cols,
row_order = chm_plot_cgc_rows,
width = 1)
chm_mut_plot_cgc <- Heatmap(chm_mut_cgc, name = "Mutation Status",
na_col = "black",
col = c("cyan3", "darkorchid"),
cluster_rows = FALSE, cluster_columns = FALSE,
show_row_names = FALSE, show_column_names = FALSE,
column_title = "Mutation Status",
column_title_gp = gpar(fontsize = 60, fontface = "bold"),
column_order = chm_plot_cgc_cols,
row_order = chm_plot_cgc_rows,
width = 1)
chm_all_plot_cgc <- chm_plot_cgc + chm_ge_plot_cgc + chm_cn_plot_cgc + chm_mut_plot_cgc
pdf("./plots_18Q3/crispr_heatmap_pearson_all_cgc.pdf", width = 100, height = 100, paper = "special")
draw(chm_all_plot_cgc, heatmap_legend_side = "right", annotation_legend_side = "right")
seekViewport("annotation_Cell_Line_SSMD")
grid.text("Cell Line SSMD", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Cas9_Activity")
grid.text("Cas9 Activity", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Culture_Type")
grid.text("Culture Type", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Primary_Tissue")
grid.text("Primary Tissue", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Primary_Disease")
grid.text("Primary Disease", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Gender")
grid.text("Gender", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
seekViewport("annotation_Source")
grid.text("Source", unit(1, "npc") + unit(2, "mm"), 0.5, gp = gpar(fontsize = 30), default.units = "npc", just = "left")
dev.off()4 Per-gene plots
4.1 Grouped plots
Plot code:
crispr_all_genes <- unique(crispr_data$Hugo_Symbol)
crispr_all_groups <- lapply(crispr_all_genes, makeCRISPRgrob)
crispr_all_groups_paths <- paste0(crispr_all_genes, "_crispr_grouped.png")
pwalk(list(crispr_all_groups_paths, crispr_all_groups), ggsave, path = "./plots_18Q3/crispr_grouped_plots", dpi = 300, width = 12, height = 12, units = "in")Data management:
crispr_all_bygene <- paste0(list.files("./plots_18Q3/crispr_grouped_plots", full.names = TRUE))
names(crispr_all_bygene) <- str_replace_all(crispr_all_bygene, c("_crispr_grouped.png" = "", "./plots_18Q3/crispr_grouped_plots/" = ""))Plots:
crispr_bygene_order <- c(intersect(crispr_signif$Hugo_Symbol, names(crispr_all_bygene)), setdiff(names(crispr_all_bygene), crispr_signif$Hugo_Symbol))
crispr_all_bygene <- crispr_all_bygene[match(crispr_bygene_order, names(crispr_all_bygene))]
bsselect(crispr_all_bygene, type = "img", selected = "KRAS", live_search = TRUE, show_tick = TRUE, height = 300, frame_height = 275)4.1.1 KRAS
makeCRISPRgrob("KRAS")
4.1.2 TP53
makeCRISPRgrob("TP53")
4.2 Sorted lineage plots
4.2.1 KRAS
makeCRISPRlinplot("KRAS")
4.2.2 TP53
makeCRISPRlinplot("TP53")
4.2.3 BRAF
makeCRISPRlinplot("BRAF")
4.2.4 NRAS
makeCRISPRlinplot("NRAS")
5 References
Barretina, J., Caponigro, G., Stransky, N., Venkatesan, K., Margolin, A. A., Kim, S., … Garraway, L. A. (2012). The Cancer Cell Line Encyclopedia enables predictive modelling of anticancer drug sensitivity. Nature, 483(7391), 603–607. https://doi.org/10.1038/nature11003
Broad Institute Cancer Dependency Map; Cancer Data Science (2018): Cancer Dependency Map, CRISPR Avana dataset 18Q3 (Avana_public_18Q3). figshare. Fileset. doi:10.6084/m9.figshare.6931364.v1
Consortium, T. C. C. L. E., & Consortium, T. G. of D. S. in C. (2015). Pharmacogenomic agreement between two cancer cell line data sets. Nature, 528(7580), 84–87. https://doi.org/10.1038/nature15736
Data Science, Cancer (2018): DEMETER2 data. figshare. Fileset. doi:10.6084/m9.figshare.6025238.v2
Doench, J. G., Fusi, N., Sullender, M., Hegde, M., Vaimberg, E. W., Donovan, K. F., … Root, D. E. (2016). Optimized sgRNA design to maximize activity and minimize off-target effects of CRISPR-Cas9. Nature Biotechnology, 34(2), 184–191. https://doi.org/10.1038/nbt.3437
Meyers, R. M., Bryan, J. G., McFarland, J. M., Weir, B. A., Sizemore, A. E., Xu, H., … Tsherniak, A. (2017). Computational correction of copy-number effect improves specificity of CRISPR-Cas9 essentiality screens in cancer cells. Nature Genetics, 49(12), 1779–1784. https://doi.org/10.1038/ng.3984
McFarland, J. M., Ho, Z. V., Kugener, G., Dempster, J. M., Montgomery, P. G., Bryan, J. G., … Tsherniak, A. (2018). Improved estimation of cancer dependencies from large-scale RNAi screens using model-based normalization and data integration. https://doi.org/10.1101/305656
6 Session information
print(sessionInfo())## R version 3.5.0 (2018-04-23)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] grid parallel stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] bindrcpp_0.2.2 circlize_0.4.4 ComplexHeatmap_1.18.1
## [4] bsselectR_0.1.0 caret_6.0-80 lattice_0.20-35
## [7] reshape2_1.4.3 devtools_1.13.6 glmnet_2.0-16
## [10] foreach_1.4.4 Matrix_1.2-14 broom_0.5.0
## [13] gridExtra_2.3 kableExtra_0.9.0 matrixStats_0.54.0
## [16] ggrepel_0.8.0 forcats_0.3.0 stringr_1.3.1
## [19] dplyr_0.7.6 purrr_0.2.5 readr_1.1.1
## [22] tidyr_0.8.1 tibble_1.4.2 tidyverse_1.2.1
## [25] plyr_1.8.4 CePa_0.6 data.table_1.11.4
## [28] rowr_1.1.3 ggpubr_0.1.7.999 magrittr_1.5
## [31] ggplot2_3.0.0 ggsignif_0.4.0 NMF_0.21.0
## [34] Biobase_2.40.0 BiocGenerics_0.26.0 cluster_2.0.7-1
## [37] rngtools_1.3.1 pkgmaker_0.27 registry_0.5
##
## loaded via a namespace (and not attached):
## [1] colorspace_1.3-2 rjson_0.2.20 class_7.3-14
## [4] rprojroot_1.3-2 GlobalOptions_0.1.0 pls_2.6-0
## [7] rstudioapi_0.7 DRR_0.0.3 prodlim_2018.04.18
## [10] lubridate_1.7.4 xml2_1.2.0 splines_3.5.0
## [13] codetools_0.2-15 doParallel_1.0.11 robustbase_0.93-1.1
## [16] knitr_1.20 RcppRoll_0.3.0 jsonlite_1.5
## [19] gridBase_0.4-7 ddalpha_1.3.4 kernlab_0.9-26
## [22] sfsmisc_1.1-2 graph_1.58.0 compiler_3.5.0
## [25] httr_1.3.1 backports_1.1.2 assertthat_0.2.0
## [28] lazyeval_0.2.1 cli_1.0.0 htmltools_0.3.6
## [31] tools_3.5.0 igraph_1.2.1 gtable_0.2.0
## [34] glue_1.3.0 Rcpp_0.12.18 cellranger_1.1.0
## [37] nlme_3.1-137 iterators_1.0.10 timeDate_3043.102
## [40] gower_0.1.2 rvest_0.3.2 DEoptimR_1.0-8
## [43] MASS_7.3-50 scales_0.5.0 ipred_0.9-6
## [46] hms_0.4.2 RColorBrewer_1.1-2 yaml_2.1.19
## [49] memoise_1.1.0 rpart_4.1-13 stringi_1.2.4
## [52] highr_0.7 bibtex_0.4.2 shape_1.4.4
## [55] lava_1.6.2 geometry_0.3-6 rlang_0.2.1
## [58] pkgconfig_2.0.1 evaluate_0.11 bindr_0.1.1
## [61] labeling_0.3 htmlwidgets_1.2 recipes_0.1.3
## [64] cowplot_0.9.3 CVST_0.2-2 tidyselect_0.2.4
## [67] R6_2.2.2 dimRed_0.1.0 pillar_1.3.0
## [70] haven_1.1.2 withr_2.1.2 nnet_7.3-12
## [73] survival_2.42-6 abind_1.4-5 modelr_0.1.2
## [76] crayon_1.3.4 rmarkdown_1.10 GetoptLong_0.1.7
## [79] readxl_1.1.0 Rgraphviz_2.24.0 ModelMetrics_1.1.0
## [82] digest_0.6.15 xtable_1.8-2 stats4_3.5.0
## [85] munsell_0.5.0 viridisLite_0.3.0 magic_1.5-8